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Abstract 

In contrast to the prevailing view in the literature, it is shown that even extremely stiff sets 
of ordinary differential equations may be solved efficiently by explicit methods if limiting 
algebraic solutions are used to stabilize the numerical integration. The stabilizing algebra 
differs essentially for systems well-removed from equilibrium and those near equilibrium. 
Explicit asymptotic and quasi-steady-state methods that are appropriate when the system is 
only weakly equilibrated are examined first. These methods are then extended to the case 
of close approach to equilibrium through a new implementation of partial equilibrium ap- 
proximations. Using stringent tests with astrophysical thermonuclear networks, evidence 
is provided that these methods can deal with the stiffest networks, even in the approach 
to equilibrium, with accuracy and integration timestepping comparable to that of implicit 
methods. Because explicit methods can execute a timestep faster and scale more favorably 
with network size than implicit algorithms, our results suggest that algebraically-stabilized 
explicit methods might enable integration of larger reaction networks coupled to fluid dy- 
namics than has been feasible previously for a variety of disciplines. 



Key words: ordinary differential equations, reaction networks, stiffness, reactive flows, 
nucleosynthesis, combustion 

PACS: 02.60.Lj, 02.30.Jr, 82.33.Vx, 47.40, 26.30.-k, 95.30.Lz, 47.70. -n, 82.20.-w, 
47.70.Pq 



Email address: guidry@utk . edu (Mike Guidry). 
Preprint submitted to Elsevier 



22 December 201 1 



1 Introduction 



Problems from many disciplines require solving large coupled reaction net- 
works. Representative examples include reaction networks in combustion chem- 
istry [1], geochemical cycling of elements [2], and thermonuclear reaction networks 
in astrophysics H3I4H . The differential equations used to model these networks usu- 
ally exhibit stiffness, which arises from multiple timescales in the problem that 
differ by many orders of magnitude H1I5I6I7II . Sufficiently complex physical sys- 
tems often involve important processes operating on widely-separated timescales, 
so realistic problems tend to be at least moderately stiff. Some, such as astrophys- 
ical thermonuclear networks, are extremely stiff, with 10-20 orders of magnitude 
between the fastest and slowest timescales in the problem. Our concern here is 
with stiffness as a numerical issue, but we remark that stiffness can have important 
physical implications because complex processes often function as they do pre- 
cisely because of the coupling of very slow and very fast scales within the same 
system. 

Books on numerical and computational methods routinely state H1I6I7H that stiff 
systems cannot be integrated efficiently using explicit finite-difference methods be- 
cause of stability issues: for an explicit algorithm, the maximum stable timestep is 
set by the fastest timescales, even if those timescales are peripheral to the main phe- 
nomena of interest. The standard resolution of the stiffness problem uses implicit 
integration, which is stable for stiff systems but entails substantial computational 
overhead because it requires the inversion of matrices at each integration step. Be- 
cause of the matrix inversions, implicit algorithms tend to scale from quadratically 
to cubically with network size unless favorable matrix structure can be exploited. 
Thus, implicit methods can be expensive for large networks. 

A simple but instructive example of stiffness is provided by the CNO cycle for 
conversion of hydrogen to helium, which powers main-sequence stars more mas- 
sive than the Sun (Fig. [[])• In the CNO cycle the fastest rates under characteristic 
stellar conditions are /3 -decays with half-lives ~ 100 seconds, but to track the com- 
plete evolution of main- sequence hydrogen burning may require integration of the 
network for hydrogen burning over timescales as large as billions of years (~ 10 16 
seconds). If one tries to implement this integration using explicit forward differ- 
encing, the largest stable integration timestep will be set by the fastest rates and 
will be of order 10 2 seconds. Thus 10 14 or more explicit integration steps could 
be required to integrate the CNO cycle to hydrogen depletion. Conversely, typical 
implicit integration schemes can take stable and accurate timesteps equal to 1-10% 
of the local time over most of the integration range, and would compute the above 
numerical integration in a few hundred implicit steps. By virtue of examples such 
as this, it is broadly accepted that explicit methods are not viable for stiff networks. 
To quote the authoritative reference Numerical Recipes Q, "For stiff problems we 
must use an implicit method if we want to avoid having tiny stepsizes." 

Our main interest lies in simulations where the reaction network is one part 
of a broader problem. Let us take as representative astrophysical thermonuclear 
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Fig. 1. The CNO (carbon-nitrogen-oxygen) cycle. On the left side the main branch of the 
cycle is illustrated with solid arrows and a side branch is illustrated with dashed arrows. On 
the right side, the main branch of the CNO cycle is illustrated in more detail. 



networks coupled to multidimensional hydrodynamics. The hydrodynamical evo- 
lution controls network conditions (temperature, density, . . . ), and the network in- 
fluences the hydrodynamics through energy production and modification of com- 
position. Solution of large networks by the usual means is costly in this context 
and even ambitious simulations use only small networks, or replace the network 
entirely by parameterization. Then a more realistic network is used in a separate 
"post-processing" step, where fixed hydrodynamical profiles computed in the orig- 
inal simulation specify the variation of temperature and density with time. Such 
approximations are especially at issue for problems like Type la supernovae, which 
are 3D asymmetric explosions powered by a complex reaction network releasing 
energy greater than that of a large galaxy on a timescale of order 1 s. 

Astrophysical reaction networks have been used to illustrate, but problems in 
various fields exhibit a similar complexity. For example, in astrochemical kinetics 
large chemical evolution networks must be modeled in dynamical environments 
such as contracting molecular clouds, or in combustion chemistry the burning net- 
works are strongly coupled to dynamical simulations of the air and fuel mixture. 
Realistic networks in all such applications may be quite large. Modeling combus- 
tion of larger hydrocarbon molecules or soot formation can require hundreds to 
thousands of reacting species with up to 10,000 reactions [1J, and realistic net- 
works for supernova explosions imply hundreds to thousands of nuclear isotopes 
with tens of thousands of reaction couplings [0. In all such cases current techniques 
do not permit the coupling of realistic reaction networks to the full dynamics and 
highly- schematic reaction networks are used in even the most realistic contempo- 
rary simulations. 



2 Stiffness under Equilibrium and Non-Equilibrium Conditions 

The pessimism engendered by the Introduction notwithstanding, it would be 
highly desirable to integrate large, complex networks by explicit means because 
explicit algorithms are simple and economical, and scale favorably with network 
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size. In principle, this could be accomplished by identifying conditions under which 
some equations in the network have an approximate analytical solution and using 
that information to remove algebraically the stiffest components from the general 
numerical solution, thereby replacing the original network with an approximation 
that permits larger stable explicit-integration timesteps. To that end, the first task 
is to understand clearly the nature of the stiffness that we wish to remove from the 
equations. 



2.1 Varieties of Stiffness 



There are several fundamentally different sources of stiffness instability in re- 
action networks that are often not clearly distinguished in the literature. Pedagogi- 
cally, discussions often emphasize an instability associated with small populations 
becoming negative because of an overly-ambitious numerical integration step, thus 
converting stable decaying-exponential terms into unstable growing exponentials. 
However, there are other instabilities that can occur, even when no population vari- 
ables become negative in an integration step. In these types of instabilities one ends 
up having to take the difference of large numbers to obtain a result very near zero. 
The numerical errors that ensue in a standard explicit approach can then accumulate 
rapidly and destabilize the network, even before any populations become negative. 
This is still a stiffness instability because the problem results ultimately from a 
numerical integration attempting to cope with highly-disparate timescales, but the 
origin of these timescales is different from that discussed above. In this case the 
disparate timescales are the fast reactions driving the system to equilibrium con- 
trasted with the slow timescale associated with equilibrium itself (which tends to 
infinity). To illustrate, consider the form of the equations that we desire to solve: 



dt 

= +#+...).- t/r +fi +■■■).- 
= (A + - fi)i + (fi -fi)i+-= Lift 



fj 



(i) 



where the yt(i = 1 . . .N) describe the dependent variables (species abundances for 
our examples), t is the independent variable (time in our examples), the fluxes be- 
tween species i and j are denoted by (/,•)/, and the sum for each variable i is over 
all variables j coupled to i by a non-zero flux (jf)i. For an N -species network 
there will be N such equations in the populations yi, with these equations generally 
coupled to each other because of the dependence of the fluxes on the different yj. 

In Eq. CD several different ways to group the terms on the right side have been 
indicated, with the first line representing a decomposition into total flux in and out 
of species i and the third line separating the total flux into in and out contributions 
from individual reactions. The utility of these alternative decompositions will be 
elaborated further below. 
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2.2 The Approach to Equilibrium 

In this discussion we shall employ the term "equilibration" broadly to mean 
that evolution of the network is being influenced strongly by nearly-canceling terms 
on the right sides of the differential equations (Q~|). Two qualitatively different equi- 
librium conditions may be distinguished that require different algebraic approaches 
to stabilization. 

2.2.1 Macroscopic Equilibration 

One class of equilibrium conditions results if F t + — F-~ — > (asymptotic) or 
F- + — Ff — > constant (steady-state). Let us refer to these conditions as a macro- 
scopic equilibration, since they are statements about entire differential equations 
in Eq. (UJ). We shall introduce asymptotic and quasi-steady-state approximations 
exploiting these conditions that remove whole differential equations from the nu- 
merical integration for a network timestep, replacing them with algebraic solutions. 
Such approximations still integrate the full original set of equations, but they reduce 
the number of equations integrated numerically by forward difference. This helps 
with stiffness because integration of some or all of the stiffest equations is replaced 
by a stable analytical solution, and any equations that remain to be integrated nu- 
merically tend to have smaller disparities in timescales and thus less stiffness. 

2.2.2 Microscopic Equilibration 

In Eq. (OQ), F} and Ff for a species i each consist of various terms depending 
on the other populations in the network. Groups of individual terms on the right 
side of Eq. (Q~|) may come into equilibrium (so that the sum of their fluxes tends 
to zero), even if the macroscopic conditions for equilibration for the entire equa- 
tion are not satisfied. We shall term this microscopic equilibration. Then we may 
consider an algebraic approximation that removes groups of such terms from the 
numerical integration by replacing their sum of fluxes with zero. This does not 
reduce the number of equations integrated numerically, but can reduce their stiff- 
ness by removing terms with fast rates, thereby reducing the disparity between the 
fastest and slowest timescales in the system. Such considerations will be the basis 
of the partial equilibrium methods to be discussed in §[5J 

3 Reaction Networks in the Context of Larger Systems 

Assume the coupling of reaction networks to a larger system to be done us- 
ing operator splitting, where the larger system (the hydrodynamical solver in our 
examples) is evolved for a timestep holding the parameters computed from the net- 
work in the previous step constant, and then the network is evolved over the time 
corresponding to the hydrodynamical timestep while holding the new variables cal- 
culated in the hydrodynamical timestep constant. This places two strong constraints 
on methods: 
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(1) The network must be advanced with new initial conditions for each hydrody- 
namical step. Thus, algorithms must initialize simply and quickly. 

(2) Integration of the network must not require time substantially larger than that 
for the corresponding hydrodynamical timestep. 

Let us elaborate further on this second point. Taking the multidimensional, adaptive- 
mesh, explicit hydrodynamical code Flash (8]| applied to Type la supernova ex- 
plosions as representative, we estimate in Ref. BUl that a simulation with realis- 
tic networks can be done in a tractable amount of time if the algorithm can take 
timesteps At over the full hydrodynamical integration range that average at least 
(0.01 — 0.001) x t, where t is the elapsed time in the integration. Timesteps of this 
size are possible with implicit and semi-implicit algorithms, but those methods are 
inefficient at computing each timestep in large networks; explicit methods can com- 
pute a timestep efficiently, but timesteps this large are unthinkable with a normal 
explicit algorithm because they are typically unstable. In this paper we shall demon- 
strate explicit integration methods that realize such large integration timesteps in a 
variety of extremely stiff examples. This alters essentially the discussion of whether 
explicit methods, with their faster computation of timesteps and more favorable 
scaling with network size, are practical for large, stiff networks. The general ap- 
proximations to be discussed have been implemented before H1I10I11I12L but we 
shall find that our implementations are much more successful than previous ap- 
plications to extremely stiff networks; accordingly we shall reach rather different 
conclusions about these methods than those of earlier publications. 



4 Algebraic Stabilization of Solutions Far from Equilibrium 

Let us address first explicit approximations that are appropriate if the system is 
far from microscopic equilibrium. (Methods to determine whether this condition is 
satisfied will be discussed in §[5J) 

4. 1 Explicit Asymptotic Approximations 

The differential equations to be solved take the form given by Eq. £[]). Gener- 
ally, F t + and Fj~ for a given species i each consist of a number of terms depending 
on the other populations in the network. For the networks that we shall consider the 
depletion flux for the population y ( - will be proportional to y ; , 

Fr = (k i l +k i 2 + ...+ki n )y l =k l y l , (2) 

where the k'j are rate parameters (in units of time -1 ) for each of the m processes 
that can deplete y/, which may depend on the populations yj and on system vari- 
ables such as temperature and density. The characteristic timescales Tj = 1 / k'j will 
vary over many orders of magnitude in the systems of interest, meaning that these 
equations are very stiff. From Eq. © the effective total depletion rate k l for y ; at a 
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given time, and a corresponding timescale z l , may be defined as 



F~ 1 

yi k l 

permitting Eq. © to be written as 

y, = Urt- d 4\- (4) 



k' \ 1 dt / 
Thus, in a finite-difference approximation at timestep t n 

Ff(tn) 1 dy t 



yi{t n ) 



k*(t n ) k'(t n ) dt 



(5) 



The asymptotic limit for the species i corresponds to the condition F t + ~ F- , im- 
plying from Eq. © that dyi/dt ~ 0. In this limit Eq. © gives a first approximation 

y>i and local error respectively, for )>;(£«), 

For small dyi/dt a correction term may be obtained by writing the derivative term 
in Eq. © as 



dt At 



(yi(t n ) -yi(t n -i)) + ^ (e^-e^ + o(At), (7) 

where 0(x) denotes terms of order x. Substitution in Eq. © then gives 

^ = F i-w^ iyM - yi '"- l) ~ ■ 

(2) (2) 

and setting y(t n ) = y n and solving for y n gives 

where a term of order (Af) 2 /(1 + k(t n )At) has been discarded. This approxima- 
tion is expected to be valid for large kAt. Another approach is to use a predictor- 
corrector scheme within such an asymptotic approximation [1]. However, it has 
been shown [|9]| that these asymptotic approximations give rather similar results 
for the networks that we shall test, so the simple formula © will be adequate for 
the present discussion. The mathematical and numerical properties of asymptotic 
approximations have been explored previously in Refs. [Ill 101 111 1211 . 

To implement an asymptotic algorithm we define a critical value K of kAt and 
at each timestep cycle through all populations and compute the product k'At for 
each species i using Eq. © and the proposed timestep At. Then, for each species i 
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(1) If k' At < K, the population is updated numerically by the explicit Euler method. 

(2) If kAt > K, the population is updated algebraically using Eq. ([8]). 
Formally explicit integration is expected to be stable if k l At < 1 and potentially 
unstable if k l At > 1 (see the discussions in Refs. [1] and [9]). Therefore, k = 1 has 
been chosen for all examples presented here. 



4.2 Quasi-Steady-State Approximations 

An alternative explicit algebraic solution to the coupled differential equations is 
possible using the quasi- steady- state (QSS) approximations developed by Mott and 
collaborators H11I12H . which drew on earlier work in Refs. H13I14I15H . Following 
Mott et al H11I12L first notice that Eq. (QQ), expressed in the form dy/dt = F + (t) — 
k(t)y(t) using Eq. © and with indices dropped for notational convenience, has the 
general solution 

y (t)=y Q e- kt + ^ T {\-e- kt ), (9) 

provided that k and F + are constant. The QSS method then uses this equation as 
the basis of a predictor-corrector algorithm in which a prediction is made using 
only initial values and then a corrector is applied that uses a combination of initial 
values and values computed using the predictor solution. In terms of a parameter 
a(r) defined by 

, , 160r 3 + 60r 2 + llr+l 
a{T) = 360r3 + 60^ + 12r+r (10) 

where r = 1 /kAt, we adopt a predictor y p and corrector y c given by H11I12I 

y y + l + a°k°At y y + l + MAf' K } 

where a is evaluated from Eq. (flOl ) with r = l/k°At, an average rate parameter is 
defined by k = | (k° + k v ), a is specified by Eq. (flOl) with r = 1 /kAt, and 

F+ = 6cF+ + (l-a)F+. 

The corrector can be iterated if desired by using y c from one iteration step as the 
y p for the next iteration step. We implement an explicit QSS algorithm based on 
the predictor-corrector (fTTT) essentially in parallel with that described above for the 
asymptotic method, except that in applying the QSS algorithm all equations are 
treated in QSS approximation, rather than dividing the equations into a set treated 
by explicit forward difference and a set treated analytically |fT6ll . 
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5 Algebraic Stabilization in the Approach to Equilibrium 

We shall find that the asymptotic and quasi-steady-state methods described in 
preceding sections work well for macroscopic equilibration, but are highly inef- 
ficient for microscopic equilibration. Thus, these methods must be augmented by 
a means to remove stiffness associated with the approach to (microscopic) equi- 
librium if they are to be applicable to a broad range of problems. In this section 
approximations to stabilize explicit integration in the presence of microscopic equi- 
libration are developed. This development draws heavily on the partial equilibrium 
work of David Mott [12] as a starting point, but we shall extend these methods and 
find much more favorable results for extremely stiff networks than those obtained 
in the pioneering work of Mott and collaborators [fTTl . 

Partial equilibrium (PE) methods examine source terms / ; + and fj~ for individ- 
ual reaction pairs in the network — not the composite fluxes F- + and F- that are the 
basis for asymptotic and QSS approximations — for approach to equilibrium. When 
a fast reaction pair nears equilibrium its source terms are removed from the direct 
numerical integration in favor of an equilibrium algebraic constraint. Reactions not 
in equilibrium still contribute to the fluxes for the numerical integrator, but once 
fast reactions are decoupled from the numerical integration the remaining system 
typically becomes (much) less stiff. Consider a representative 2-body reaction and 
its source term f ab ^ cd , 

a+b^c+d fabled = ± (k&ayb - k ry c yd ), (12) 

where the y, denote population variables for the species i and the ks are rate pa- 
rameters for forward (kf) and reverse (k T ) reactions. Considered in isolation, the 
reaction pair of Eq. (fT2l may be deemed to be in equilibrium if fabled = 0- It is 
useful to introduce the idea of partial equilibrium (PE), where at a given time some 
reaction pairs have / = and some have / 7^ 0. The evolution of the system is then 
determined primarily by those reactions for which / 7^ 0, but since the system is 
coupled the / 7^ reactions will perturb the f = reaction pairs so that for those 
reactions near equilibrium / = — > f ~ 0. Let us now introduce a set of definitions 
and concepts that will allow us both to quantify and use to our advantage the partial 
equilibrium condition / ~ 0. 

5. 1 Conserved Scalars and Progress Variables 

The reaction pair a + b c + d appears at first glance to have four characteristic 
timescales associated with the rate of change for the four populations y a , yb, y c , and 
yd, respectively. However, it is clear that the following three constraints apply to 
this reaction 

y a -yb = c\ y a +y c = c2 y a +yd = c3, (13) 
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where the constants c/ may be evaluated by substituting the initial abundances into 
Eq.O, 

ci=y° a -y° b c 2 =y°+y Q c <* =35 + 35- 
Losing one a in the reaction a + b ^ c + d requires the simultaneous loss of one 
b, so their difference must be constant, as implied by the first of Eqs. (PT3l) . and 
every loss of one a produces one c and one d, implying the second and third of Eqs. 
(TT3T) . The left sides in the equations (fT3l) are examples of conserved scalars |fT2|. 
which are constant by virtue of the structure of the equations, not by any particular 
dynamical assumptions. The differential equation describing the evolution of y a is 

dy a 2 

—r = -hy a yb + k r y c yd = ay a +by a +c 1 (14) 

where Eq. (fT3l) has been used and 

a = k r — kf b = —k r (c2 + C3) + kfC\ c = k^c^,. 

We shall demonstrate below that the approach to equilibrium for any 2-body reac- 
tion pair can be described, and the approach to equilibrium for any 3-body reaction 
pair can be approximated, by a differential equation of this form. In terms of the 
quantity 

q = 4ac-b 2 : (15) 

the solutions to Eq. (fl4l) of interest in the present context correspond to a 7^ and 
q < 0, and take the form H12I18H 

M 1 (ua I — l + 0exp(- A /=gr) > \ _ 2ay + b + y/=q~ 

ya{t) = -—[b + yj-q- 7=-r = ^ 1=- (16) 

2a \ 1 — ^)exp(— y/— qt) J layo + b — ^—q 



The equilibrium solution then corresponds to the limit t — > °° of Eq. (1161) . 

ya = y^ = -^(b + x/=q). (17) 

Once y a has been determined the constraints (fT3l) may be used to determine the 
other abundances. For example, at equilibrium 

y~b{t)=y a {t)-c\ y c (t) = c2-y a (t) ydif) = <%-%{$). 

It is often convenient to introduce a variable A that is the difference between the 
values of the v, at the beginning of the timestep and their current values 

y a = y°a-A. yb = y° b -h y c =y° c + ^ y d = y d +^ (18) 

The new variable A is termed a progress variable for the reaction characterized by 
fabled and satisfies 

= fabled A = X (t = 0) = 0. (19) 
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Fig. 2. (a) Time evolution of the solution ( fT6l ) assuming constant a, b, and c. The charac- 
teristic timescale for approach to equilibrium (f2Ql > is labeled x and the equilibrium value 
of y(t) defined by Eq. ( fTTl ) is denoted by y. To illustrate we have assumed the initial value 
yo = 0. For times considerably larger than x the general solution (fT6l ) saturates at the equi- 
librium solution (fTTT ). (b) Behavior of the exponential factor in Eq. ( fTBT ). 



Thus the approach of a + b ^ c + d to equilibrium is controlled by a single dif- 
ferential equation (lT4l) that can be expressed in terms of either a single one of the 
abundances yi, or the progress variable A. The general solution of this equation is 
of the form given by Eq. (fT6l) . with the time dependence residing dominantly in 
the exponentials. Therefore, the rate at which a + b # c + d evolves toward the 
equilibrium solution (fTTT) is governed by a single timescale 




which is illustrated in Fig. [2l Whether a reaction pair is near equilibrium at time t 
may then be determined by requiring that 

< £l (21) 

for each of the species i that is involved in the reaction pair, where yi(t) is the 
actual abundance, yi is the equilibrium abundance determined from Eq. (fTTT) . and 
£, is a user-specified tolerance that could depend on / but will be taken to be the 
same for all species in the examples discussed to be here. Alternatively, the equi- 
librium timescale (1201) may be compared with the current numerical timestep to 
determine whether a reaction is near equilibrium: if x is much smaller than the 
timestep, we may expect that equilibrium can be established and maintained in suc- 
cessive timesteps, even if it is being continually disturbed by other non-equilibrated 
processes. 
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5.2 Reaction Vectors 



A partial equilibrium approximation in a large network could require that thou- 
sands of reactions be examined for their equilibrium status at each timestep. Let us 
introduce a formalism, adapted from the work of Mott [fT2l . that allows examina- 
tion of partial equilibrium criteria in a particularly efficient way by exploiting the 
analogy of a reaction network to a linear vector space. This will have two large 
advantages: (1) It will provide us with some well-established mathematical tools. 
(2) The abstraction of reaction network as linear vector space permits formulation 
of a partial equilibrium algorithm that is not strongly tied to the details of a partic- 
ular problem, thus aiding portability within and across disciplines. 

We begin by expressing the concentration variables for the n species A; in a 
network as components of a composition vector 



which lies in an n-dimensional vector space 4>. The components v ; are proportional 
to number densities for the species labeled by Au so a specific vector in this space 
defines a particular composition. Any reaction in the network can then be written 
in the form YJ!=\ a ; A ; ^ YH=\ bi^i for some sets of coefficients {a,} and {£>,}. The 
coefficients on the two sides of a reaction may be used to define a vector r E 4> with 
components 

r=(bi-ai,b2-a2,...,b n -a n ), (23) 

that specifies how the composition may change because of the reaction. For exam- 
ple, consider the main part of the CNO cycle illustrated on the right side of Fig. [TJ 
Choosing a basis {/?, a, 12 C, 13 N, 13 C, 14 N, 15 0, 15 N}, the reaction p + 15 N -»■ a + 
12 C then has a reaction vector r with components (— 1 , 1 , 1 , 0, 0, 0, 0, — 1 ) . 

For a network with three or fewer species the corresponding linear vector space 
can be displayed geometrically. For example, Fig.[3]illustrates a network containing 
the isotopes { 4 He, 8 Be, 12 C} and four reaction vectors corresponding to the reaction 
pairs 2a # 8 Be and a + 8 Be # 12 C. Larger networks cannot be visualized so easily, 
but their algebraic properties remain completely analogous to those of a simple 
network like that illustrated in Fig. [3] 

5.3 Conservation Laws 

Given an initial composition y = (y^y®, ■ ■ -Jn)' a single P au " of reactions la- 
beled by i can produce a composition y = y + so for a set of k possible 
reactions 



y=(y\,y2,y3,---yn), 



(22) 



k 



(24) 
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Fig. 3. (a) Four reaction vectors for a simple 3-species network corresponding to the astro- 
physical triple-a process that converts 4 He to 12 C in red giant stars. 

Define a time-independent vector cG<J> that is orthogonal to each of the k vectors 
r t in Eq. (1241) . Then from Eq. (|24l) . c-y = c-y and therefore 

n n 

= H Ci y'o = c ° nstant 5 ( 25 ) 

!=1 ;=1 

so any vector c orthogonal to the reaction vectors r\, r^, ■■■T\ gives a linear com- 
bination of abundances that is invariant under all actions of the network. Equation 
(|25l) is a conservation law of the system that follows purely from the structure of the 
network and thus is valid irrespective of dynamical conditions in the network. The 
conserved quantities may be determined by solving the equation r ■ c = 0, where r 
is a matrix with rows composed of all reaction vectors (|23T) . for the vectors c^. 

Thus, once a basis has been chosen for a specific reaction network it is straight- 
forward to enumerate all of its possible compositions in terms of a set of vectors y h 
all possible reactions in terms of a set of vectors ry, and the conservation laws that 
follow from the structure of the network in terms of a set of vectors c^. This abstract 
description of a reaction network retains a reference to a specific physical problem 
only through the choice of basis and the choice of allowed reaction vectors. 

5.4 Reaction Group Classes 

It will prove useful to associate inverse reaction pairs in reaction group classes 
(reaction groups, or RG for short), in which all reactions of a group share the same 
reaction vector r, up to a sign. To illustrate we employ the reaction classifications 
used in the REACLIB library |[T9ll that are illustrated in Table [TJ In this classi- 
fication the reactions of importance in nuclear astrophysics are assigned to eight 
categories, depending on the number of nuclear species on the two sides of the 
reaction equation. From this classification we see that there are five independent 
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Table 1 

Reaction classes in the REACLIB lfl9ll library 



Place 


rvCaCLlUll 


Tioc/m*i r^l ii"vn / ~\ t- pvQinnlp 
L-/C»CI1UL1U11 UI CAdlllUlC 


1 
1 


ct — )" b 


p-oecay or e capture 





a r U i C 


r IlUlUUlMllieglallUll t£ 


J 


o v K i f i rl 


12 r — v ^ry 


4 


a + b — >■ c 


Capture reactions 


5 


a + b -> c + d 


Exchange reactions 


6 


a + b^c + d + e 


2 H + 7 Be^ 1 H + 2 4 He 


7 


a + b— ^c + d + e + f 


3 He + 7 Be^2 1 H + 2 4 He 


8 


a + b + c^d (+ e) 


Effective 3-body reactions 



Table 2 



Reaction group classes 



Class 


Reaction pair 


REACLIB class pairing 


A 


a^b 


1 with 1 


B 


a + b # c 


2 with 4 


C 


a +b + c ^ d 


3 with part of 8 


D 


a + b c + d 


5 with 5 


E 


a + b^c + d + e 


6 with part of 8 



ways to combine the reactions of Tabled] into reversible reaction pairs, leading to 
the reaction group classification illustrated in Table [2] For example, reaction group 
B consists of reactions from REACLIB reaction class 2 (a — > b + c) paired with 
their inverse reactions (b + c — > a), which belong to REACLIB reaction class 4. 

For each reaction group class the characteristic differential equation governing 
the reaction pair considered in isolation is of the form given by Eq. (fl4"l) . dy/dt = 
ay 2 + by + c, where y is either an abundance variable (proportional to a number 
density) for one of the reaction species, or a progress variable measuring the change 
from initial abundances associated with the reaction pair, and the coefficients a, b, 
and c are known parameters depending on the reaction rates that will be assumed 
constant within a single network timestep. An exception occurs for reaction group 
classes C and E, where there are a 3-body reactions and the most general form of the 
differential equation involves cubic terms, dy/dt = ay 3 + f5y 2 + yy + e. These "3- 
body" reactions in astrophysics are typically actually sequential 2-body reactions 
and we assume that in any timestep y(t) 3 ~ y^y{t) 2 , where is the (constant) 
value of y(t) at the beginning of the timestep. This reduces the cubic equation to 
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an effective quadratic equation of the form (fl?)) . with a = ccy^ + j8, b = y, and 
c = e. This has been found to be a very good approximation in typical astrophysical 
environments and we apply it to all 3-body reactions in the examples discussed 
here. 

5.5 Equilibrium Constraints 

If a reaction pair from a specific reaction group class is near equilibrium, there 
will be a corresponding equilibrium constraint that takes the general form ffTTTl 

n#~ ai) =^ (26) 

z=l 

where K is a ratio of rate parameters. For example, consider the reaction group class 
E pair a + b^ c + d + e in isolation, with differential equations for the populations 
y,- written in the form 

% = % = ~% = -yd = -% = -hy a yb + k r y c y d y e . 

Then at equilibrium, requiring that the forward flux — kfy a yb an d reverse flux k r y c ydy e 
in the reaction pair sum to zero gives the constraint 

y a yb _^ =K 
y c yaye h ~ 

which is of the form (1261) . 

5.6 Systematic Classification of Reaction Group Properties 

Applying the principles discussed in the preceding paragraphs to the reaction 
group classes in Table [2] gives the results summarized for reaction group classes 
A-E in Appendix lAl This reaction group classification has been developed assum- 
ing astrophysical thermonuclear networks and a particular parameterization (REA- 
CLIB) of the corresponding reaction rates. However, the illustrated methodology 
is of wider significance. First, since any reaction compilation in astrophysics could 
be reparameterized in the REACLIB format, this classification scheme provides a 
general partial-equilibrium bookkeeping that could be applied to any problem in 
astrophysics. Second, for any large reaction network in any field we may apply 
the classification techniques illustrated here to group all reactions of the network 
into reaction group classes and deduce for each reaction group class analytical ex- 
pressions for all quantities necessary for applying a PE approximation. All that is 
required is to cast the network as a linear algebra problem by choosing a set of 
basis vectors corresponding to the species of the network, and then to define the 
corresponding reaction vectors of physical importance within this space. Once that 
is done, the formalism developed here may be applied systematically. In principle 
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this classification need be developed only once for the networks of importance in 
any particular discipline, as we have illustrated here specifically for astrophysics. 

5. 7 General Methods for Partial Equilibrium Calculations 

We now have a set of tools to implement partial equilibrium approximations, 
but there are some practical issues to resolve before it is possible to make realis- 
tic calculations. To address those issues, let us now outline a specific application 
approach. Although our remarks will be illustrated by examples using astrophysi- 
cal thermonuclear networks, the methods discussed should be relevant for a much 
broader range of problems. 

5. 7. 1 Overview of Approach 

We employ the partial equilibrium method in conjunction with the asymptotic 
approximation. (There may be advantages in using quasi- steady- state methods in- 
stead of asymptotic ones, but we shall deal with that in future work.) Once the 
reactions of the network are classified into reaction groups, the algorithm has three 
steps: 

(1) A numerical integration step begins with the full network of differential equa- 
tions, but in computing fluxes all terms involving reaction groups judged to be 
equilibrated (based on criteria determined by species populations at the end of 
the previous timestep) are assumed to sum identically to zero net flux and are 
omitted from the flux summations. 

(2) A timestep At is chosen and used in conjunction with the fluxes to determine 
which species satisfy the asymptotic condition. For those species that are not 
asymptotic, the change in abundance for the timestep is then computed by or- 
dinary forward (explicit) finite difference, but for those species that are judged 
to be asymptotic the new abundance for the timestep is instead computed an- 
alytically using the asymptotic approximation. 

(3) For all species in reaction groups judged to be equilibrated at the beginning of 
the timestep, it is assumed that reactions not in equilibrium will have driven 
these populations slightly away from equilibrium during the timestep. These 
populations are then adjusted, subject to the system's conservation laws, to 
restore their equilibrium values at the end of the timestep. 

Hence the partial equilibrium part of the approximation does not reduce the num- 
ber of differential equations integrated numerically within a timestep, but instead 
removes systematically the stiffest parts of their fluxes. In contrast, the asymptotic 
approximation effectively reduces the number of differential equations integrated 
numerically by replacing the numerical forward difference with an analytic asymp- 
totic abundance for those isotopes satisfying the asymptotic condition. 

The partial equilibrium and asymptotic approaches are complementary because 
partial equilibrium can operate microscopically to make the differential equation 
for a given isotope less stiff, even if that isotope does not satisfy the asymptotic con- 
dition, while the asymptotic condition removes entire differential equations from 
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the numerical update and thus can operate macroscopically to remove stiff reaction 
components, even if they do not satisfy partial equilibrium conditions. For brevity 
we shall often refer simply to the partial equilibrium (PE) approximation, but this 
always means the partial equilibrium approximation used in conjunction with the 
asymptotic approximation, in the manner just described. 

5. 7.2 Complications in Realistic Networks 

The complication for the basic approach outlined in §15.7.11 when applied to a 
realistic network is that there may be more than one computed equilibrium value 
for a given species. This is because the equilibrium abundance from Eq. (TT71) is 
to be computed separately for each reaction group, and a given isotope often will 
be found in more than one such group. In the general case these computed equi- 
librium values can be different, since equilibrium within each reaction group is 
specified only to the tolerance implied by £, in Eq. (I2TT) . However, the equilibrium 
abundances of a given species computed for each equilibrated reaction group that 
it is a member of cannot differ substantially among themselves, since this would 
contradict the evidence supplied by the network that the reaction groups are in 
equilibrium. For a species i there is only one actual abundance y ; in the network at 
a given time, which must satisfy simultaneously the equilibrium conditions for all 
reaction groups determined to be in equilibrium, within the tolerances of Eq. (I2TT) . 

Therefore, restoration of equilibrium at the end of a numerical integration 
timestep will correspond to setting the abundance of each isotope to a compromise 
choice among each of the (similar) predicted equilibrium values for all equilibrated 
reaction groups in which it participates. This is a self-consistent approximation if 
the variation in possible equilibrium abundances remains less than the tolerances 
used to impose equilibrium in Eq. (121) . For the networks that we have tested, this 
seems generally to be fulfilled for appropriate choices of 

5. 7. 3 Specific Methods for Restoring Equilibrium 

We have investigated three methods for restoring equilibrium at the end of the 
numerical timestep: 

(1) Reimpose equilibrium ratios (Eq. (|26l) ) by Newton-Raphson iteration. 

(2) Reimpose equilibrium abundances (Eq. (flTl) ) by Newton-Raphson iteration. 

(3) Reimpose equilibrium abundances (Eq. (fTTT) ) by averaging over progress vari- 
ables. 

All three methods are described in Ref. ifTTl . where we show that for the examples 
investigated to date the third method gives essentially the same results as the first 
two and is considerably simpler to implement, since it involves no matrices and no 
iteration. Thus we outline only this method and refer the reader to Ref. IfPTll for the 
technical details on it and the other two methods. 

The basic idea is that in partial equilibrium the isotopic abundances in a reac- 
tion group evolve according to a single timescale given by Eq. (1201) . as discussed in 
§15.11 Thus, within a single reaction group, the equilibrium abundance of any one 



17 



isotope, or the progress variable for the reaction group, determines the equilibrium 
abundance of all species in the group. Furthermore, within a single reaction group 
the evolution of all species to equilibrium conserves particle number, by virtue of 
constraints such as those of Eq. (IT3T ) that are summarized for all five reaction group 
classes in Appendix lAl Let us exploit this using the progress variable A; from each 
reaction group. If all the reaction groups were independent, then equilibrium could 
be restored by requiring for all reaction groups in equilibrium A; — A,; = 0, where 
Xi denotes the equilibrium value of A, computed from Eq. (fTTT ) and relations like 
Eq. (TTSl . Once the equilibrium value of A; is computed, the equilibrium values for 
all other isotopes in the group than follow from constraints like Eq. (TTSl that are 
tabulated for all reaction group classes in Appendix lAl 

The simple considerations of the preceding paragraph are insufficient because 
the reaction groups are generally not independent: the individual network species 
may appear in more than one reaction group, as discussed in §15.7.21 The equi- 
librium condition implies that we must have approximately equal computed equi- 
librium values for an isotope participating in more than one RG, but exact equal- 
ity does not hold because of the finite tolerance £, used to test for equilibrium in 
Eq. (|2T|) . Thus, we restore equilibrium for each isotope participating in partial equi- 
librium at the end of a timestep by replacing its numerically-computed abundance 
with its equilibrium value (fT7l) . averaged arithmetically over all equilibrated reac- 
tion groups in which it participates. 

There is one final issue: evolution to equilibrium for individual reaction groups 
considered in isolation conserves particle number, but the averaging procedure in- 
troduces a small fluctuation since the average will generally differ from the indi- 
vidual yi that were computed conserving particle number. Thus, for each timestep, 
after equilibrium has been restored, we rescale all v, by a multiplicative factor that 
restores the total nucleon number to its value at the start of the timestep. 



6 A Simple Adaptive Timestepper for Explicit Integration 

For testing the asymptotic and QSS methods a simple adaptive timestepper has 
been used that is described in more detail in Ref. [0: 

(1) Compute a trial integration timestep based on limiting the change in popula- 
tions that would result from that timestep to a specified tolerance. Choose the 
minimum of this trial timestep and the timestep taken in the previous integra- 
tion step and update all populations by the asymptotic or QSS algorithms. 

(2) Check for conservation of particle number within a specified tolerance range. 
If not satisfied, increase or decrease the timestep as appropriate by a small 
factor and repeat the timestep using the original fluxes. 

Our asymptotic plus PE calculations use a modified form of this timestepper that 
de-emphasizes the second step, since the PE algorithm itself implements approxi- 
mate probability conservation. This prescription is far from optimized, but it gives 
stable and accurate results for the varied astrophysical networks that have been 
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Table 3 



Speedup factors for explicit vs. implicit timesteps [ 18 ] 



Network 


Isotopes 


Speedup F 


PP 


7 


~ 1.5 


Alpha 


16 


3 


Nova 


134 


7 


150-isotope 


150 


7.5 


365-isotope 


365 


~20 



tested. Thus it is adequate for our primary task here, which is to establish whether 
explicit methods can even compete with implicit methods for stiff networks. 

7 Comparisons of Explicit and Implicit Integration Speeds 

In the following examples explicit and implicit integration speeds will be com- 
pared using algorithms at different stages of development and optimization. Thus, 
the codes cannot yet be simply compared head to head. We shall make a simple 
assumption that for codes at similar levels of optimization the primary difference 
between explicit and implicit methods would be in the extra time spent in matrix 
operations for the implicit method. Thus, if the fraction of computing time spent in 
the linear algebra is / for an implicit code, it will be assumed that an explicit code 
at a similar level of optimization could compute a timestep faster by a factor of 
F = 1/(1—/). We shall then assume that for implicit and explicit codes at similar 
levels of optimization the ratio of speeds for a given problem would be F times the 
ratio of integration steps that are required by the two codes. As discussed later, this 
may underestimate the speed of a fully-optimized explicit code relative to implicit 
codes, but establishes a lower limit on how fast the explicit calculation can be. 

Estimated factors F are shown in Table [3] for the networks to be discussed, 
based on data obtained by Feger |[T8l using the implicit, backward-Euler code Xnet 
11201 with both dense solvers (LAPACK (ZD) and sparse solvers (MA28 B2l and 
PARDISO 112310 . assuming the optimal solver to be used by the implicit code for a 
given network. We see that for very small networks implicit and explicit methods 
require similar times to compute a timestep, but for larger networks the explicit 
computation of a timestep can be faster than that for the implicit code by a factor 
of 10 or more for the networks examined here. 

8 Tests of Asymptotic and QSS Algorithms 

This section presents some representative calculations for astrophysical ther- 
monuclear networks using the explicit asymptotic (Asy) and quasi- steady- state 
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Fig. 4. Integration of the pp-chains at constant temperature = 0.016 (where Tg denotes 
temperature in units of 10 9 K) and constant density 160 gem -3 , assuming solar initial 
abundances. Reaction rates were taken from the REACLIB library 1 19 ]. (a) Mass fractions 
for the asymptotic method, the QSS method, and for the standard implicit code Xnet GDI , 
(b) Integration timesteps. 

(QSS) methods. For astrophysical networks we shall replace the generic y t of Eq. CD 
with population variables common for astrophysics: the mass fraction Xj or the (mo- 
lar) abundance Yu with 



njAj 
pN A 



Y . = Xi H 
1 ~ Ai 



pN A 



(27) 



where N A is Avogadro's number, p is the total mass density, A ; is the atomic mass 
number, n ; the number density for the species i, and conservation of nucleon num- 
ber requires = 1. 



8. 1 Explicit Asymptotic and QSS Integration of the pp-Chains 



The pp-chains that power the Sun provide a spectacular illustration of stiffness 
in a simple yet physically-important network. Figure 0] exhibits integration of the 
pp-chains at a constant temperature and density characteristic of the current core of 
the Sun using the asymptotic method, the QSS method, and the implicit backward- 
Euler code Xnet [20J. We see that the asymptotic and QSS integrations give results 
for the mass fractions in rather good agreement with the implicit code over 20 or- 
ders of magnitude, and generally take timesteps dt ~ O.lt that are comparable to 
those for the implicit code over the entire range of integration. (The asymptotic 
method required 333 total integration steps, the QSS method required 286 steps, 
and the implicit method required 176 steps for this example.) The maximum sta- 
ble timestep for a standard explicit integration method, which typically may be 
approximated by the inverse of the fastest rate in the system, is illustrated by the 
dashed blue curve in Fig.Htb). Thus, at the end of the calculation the explicit inte- 
gration timesteps are ~ 10 21 times larger than would be stable in a normal explicit 
integration. The calculation illustrated in Fig. |4] takes a fraction of a second on a 
3 GHz processor with the asymptotic, QSS, or implicit methods. In contrast, from 
Fig.Stb) we estimate that a standard explicit method would require ~ 10 21 s of pro- 
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Fig. 5. (a) Mass fractions for a network under nova conditions, corresponding to the hy- 
drodynamical profile shown in (b). The calculation used the QSS method and a network 
containing 134 isotopes coupled by 1531 reactions, with rates taken from the REACLIB 
library |fl9l and initial abundances enriched in heavy elements [24]. 

cessor time to compute the pp-chains to hydrogen depletion, which is a thousand 
times longer than the age of the Universe. 

This example is a simple one, but the results call into question most of what has 
been said in the literature concerning the use of explicit methods for stiff systems. 
This network is about as stiff as one will find in any scientific application, with a 
maximum integration timestep that is 21 orders of magnitude larger than the inverse 
of the fastest rate in the system. Yet both the explicit asymptotic and explicit quasi- 
steady- state methods have integrated it with an efficiency and accuracy comparable 
to that of a standard implicit algorithm. 

8.2 Nova Explosions 

The network of the preceding example was extremely stiff, but small. Let us 
now examine a case involving an extremely stiff but much larger network. Figure 
(3a) illustrates a calculation using the explicit QSS algorithm with a hydrodynami- 
cal profile 11251 displayed in Fig.[5tb) that is characteristic of a nova outburst. Since 
there are so many mass-fraction curves we do not attempt to display a direct com- 
parison with an asymptotic or implicit calculation, but note that the agreement is 
rather good, and the total integrated energy release corresponding to the simula- 
tion of Fig. [5] was within 1% of that found for the same network using the explicit 
asymptotic approximation. 

The integration timesteps for the calculation in Fig.Oa) are displayed in Fig.(6fa). 
Once burning commences, the QSS and asymptotic solvers take timesteps that are 
from 10 6 to 10 10 times larger than would be stable for a normal explicit integration. 
The explicit QSS timesteps illustrated in Fig. [(J a) are somewhat larger than those 
of the asymptotic solver, and comparable to or greater than those for a typical im- 
plicit code: in this calculation the implicit method required 1332 integration steps, 
the asymptotic calculation required 935 steps, and the QSS method required 777 
steps. 
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Fig. 6. (a) Timesteps for integration of Fig. |5j (b) Fraction of reactions that reach partial 
equilibrium in the QSS calculation of Fig. [5] 

Given that for a network with 134 isotopes the explicit codes should be able to 
calculate an integration timestep about 7 times faster than the implicit code because 
they avoid the manipulation of large matrices (Table [3]), these results suggest that 
the explicit QSS algorithm is capable of calculating the nova network about 12 
times faster and the explicit asymptotic algorithm about 10 times faster than a state- 
of-the-art implicit code. This impressive integration speed for both the QSS and 
asymptotic methods applied to a large, extremely stiff network is possible because 
few reactions reach equilibrium during the simulation, as determined from Eq. (T2TI) 
and illustrated in Fig.[6tb). 

The ability of explicit methods to deal effectively with large networks under 
nova conditions is supported further in work reported by Feger, et al H18I26L There 
results similar the present ones were demonstrated with asymptotic methods for a 
nova simulation using a different network, different hydrodynamical profile, and 
different reaction library than those employed here. 

8.3 Tidal Supernova Simulation 

A comparison of asymptotic and QSS mass fractions and timesteps is shown in 
Fig- Eta) for an alpha network with a hydrodynamical profile illustrated in Fig.[7£b) 
that is characteristic of a supernova induced by tidal interactions in a white dwarf 
[1271 . The integration timestepping is compared for QSS, an asymptotic calculation, 
and the implicit code Xnet in Fig. Oa). We see that the timestepping for the QSS 
calculation is somewhat better than for the asymptotic code and considerably better 
than for the implicit code (242 total integration steps for the QSS calculation, 480 
steps for the asymptotic calculation, and 1185 steps for the implicit calculation). 
Since from Table [3] an optimized explicit code could compute timesteps about 3 
times as fast as an implicit code for an alpha network, we may estimate that the 
QSS code is capable of calculating this network perhaps 15 times faster, and the 
asymptotic code perhaps 7 times faster, than the implicit code. The relatively good 
timestepping for the QSS and asymptotic methods in this case again is primarily 
because almost no reactions in the network come into equilibrium over the full 
range of the calculation, as illustrated in Fig. Ob). 
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Fig. 7. Comparison of asymptotic and QSS approximation calculations for an alpha net- 
work calculation under tidal supernova conditions. Initial composition was pure 4 He and 
REACLIB rates lfl9l were used, (a) Mass fractions, (b) Hydrodynamical profile [27]. 




Fig. 8. (a) Integration timesteps and maximum stable purely-explicit timestep ~ 1 /R max for 
the calculation in Fig. [7] (b) Fraction of isotopes that become asymptotic and fraction of 
reactions equilibrated in the network for the calculation in Fig. |7] 



A calculation for the hydrodynamical profile illustrated in Fig.[7£b) for a 365- 
isotope network using asymptotic and implicit methods is illustrated in Fig. [9] The 
implicit calculation required 1455 integration steps, compared with 5778 steps for 
the asymptotic calculation. But Table [3] indicates that for this 365-isotope network 
the explicit code can calculate each timestep about 20 times faster than the im- 
plicit code, so an optimized asymptotic code should be capable of performing the 
integration in Fig. |9] perhaps 5 times faster than a state-of-the-art implicit code. 

This ability of explicit methods to deal effectively with large networks under 
tidal supernova conditions is supported further by results from Refs. H18I26L Al- 
though different networks and different reaction network rates were used in these 
references, the explicit asymptotic method was again found to be highly competi- 
tive with standard implicit methods for the tidal supernova problem. 



8.4 Non-Competitive Explicit Timesteps in the Approach to Equilibrium 



The examples shown to this point have deliberately emphasized examples from 
networks in which few reactions have become microscopically equilibrated. For 
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Fig. 9. (a) Mass fractions calculated in asymptotic approximation using asymptotic methods 
for a 365-element network with 4325 reaction couplings under tidal supernova conditions, 
corresponding to the hydrodynamical profile shown in Fig. |7Jb). Rates were taken from 
REACLIB [ 19] and the initial composition was pure 4 He. (b) The corresponding integration 
timesteps and maximum stable explicit timestep ~ 1 /7? max - 




Log time (s) Log time (s) 

Fig. 10. Comparison of asymptotic and QSS approximations for an alpha network with 
constant temperature Tg = 5 and constant density of 10 8 gem -3 , using REACLIB rates 
|[T9l and initial equal mass fractions of 12 C and 16 0. (a) Isotopic mass fractions, (b) Inte- 
gration timesteps (left axis) and fraction of reaction in partial equilibrium (right axis). The 
gray shaded area represents roughly the region that the explicit timestep profile must lie 
in to have a chance to compete with implicit methods. The different asymptotic methods 
are labeled Asy and are described in Ref. |9). The dashed green line (PE) represents the 
timestepping using the partial equilibrium methods that will be discussed in $9l 



such cases the estimated integration speed for optimized quasi-steady-state and 
asymptotic explicit methods is often comparable to, and in some cases may ex- 
ceed, that for current implicit codes. Let us now turn to an example representative 
of a whole class of networks where this is decidedly not true. The calculation in 
Fig. [10] compares various asymptotic approximations and a QSS calculation with 
an implicit calculation for an alpha-particle network (see Table© at a constant tem- 
perature and density that might be found in a deflagration-burning zone for a Type 
la supernova simulation. Two important conclusions follow from these results. 
(1) All of the QSS and asymptotic cases shown are similar, with integrated final 

energies (not shown) that lie within 1% of each other and total integration 

times that lie within 25% of each other. 
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(2) The asymptotic and QSS methods all give timesteps that are potentially com- 
petitive with implicit methods (they lie in the shaded gray area) at earlier 
times, but fall far behind for times greater than about 10 5 s. 
The reason for conclusion (2) may be seen from the solid red curve on the right 
of Fig. [TOf b). which represents the fraction of reactions that satisfy partial equilib- 
rium conditions. Asymptotic and quasi-steady-state approximations work very well 
as long as the network is well-removed from equilibrium, but as soon as significant 
numbers of reactions become equilibrated the asymptotic and QSS timestepping 
begins to lag by large margins. As we now show, partial equilibrium methods must 
be used in conjunction with asymptotic or QSS methods to recover competitive 
timestepping in the approach to equilibrium. A preview of those results is displayed 
in Fig. [TOTb). The dashed green curve labeled PE corresponds to an asymptotic 
plus partial equilibrium approximation that exhibits highly-competitive timestep- 
ping relative to that of the implicit calculation, even as the network approaches 
equilibrium. 

9 Tests of Partial Equilibrium for Some Thermonuclear Alpha Networks 

The partial equilibrium algorithm of ^5] has been tested in a variety of ther- 
monuclear alpha networks. This section gives some representative examples of 
those calculations. Because they are among the toughest numerical problems for re- 
action networks, we shall concentrate on conditions expected in Type la supernova 
explosions (temperatures from 10 7 -10 10 K and densities from 10 7 -10 9 gcm~ 3 ). 
Since the Type la explosion is triggered by a thermonuclear runaway in the degen- 
erate matter of a white dwarf, the reaction network must deal with temperatures that 
can change at rates in excess of 10 17 K/s. Such conditions lead to rapid equilibra- 
tion in systems of extremely high stiffness and represent a demanding test of any 
numerical integration method. Although our longer-term goal is application of the 
present methods to larger networks with hundreds or even thousands of isotopes, 
an alpha network provides a fairly realistic, highly-stiff system for initial tests that 
has the advantage of being small enough to provide transparency in how the algo- 
rithm functions. It should also be noted that the most ambitious calculations to date 
for astrophysical thermonuclear networks coupled to hydrodynamical simulations 
have employed alpha networks (or even more schematic ones). Reactions and cor- 
responding reaction groups are displayed in TablelU REACLIB [19] was used for 
all rates except that inverse rates for reaction groups 3, 5, and 6 are not included 
in the standard REACLIB tabulation and were taken from Ref. ll28Tl . For all calcu- 
lations the equilibrium criteria were imposed using Eq. (T2TI) . with a constant value 
£, = 0.01 for the tolerance parameter. 

9. 1 Example at Constant Temperature and Density 

A calculation for constant Tg = 5 and density of 1 x 10 7 g cm~ 3 in an alpha net- 
work is shown in Fig. [CD The calculated asymptotic plus partial equilibrium (Asy + 
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Table 4 

Reactions of the alpha network for partial equilibrium calculations. The reverse reactions 
such as 20 Ne — > a + 16 are photodisintegration reactions, 7 + 20 Ne — > a + 16 0, with the 
photon 7 suppressed in the notation since we track only nuclear species in the network. 



Group 


Class 


Reactions 


Members 


1 


C 


3a # 12 C 


4 


2 


B 


a + 12 C^ 16 


4 


3 


D 


12 C + 12 C^a + 20 Ne 


2 


4 


B 


a + 16 # 20 Ne 


4 


5 


D 


12 C + 16 0^a + 24 Mg 


2 


6 


D 


16 + 16 0^a + 28 Si 


2 


7 


B 


a + 20 Ne^ 24 Mg 


4 


8 


D 


12 c + 20 Ne _^ a + 28 S j 


2 


9 


B 


a + 24 Mg ^ 28 Si 


4 


10 


B 


a + 28 Si^ 32 S 


2 


11 


B 
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PE) mass fractions are compared with those of an implicit code in Fig.QjJa). There 
are small discrepancies in localized regions, especially for some of the weaker pop- 
ulations, but overall agreement is rather good. Mass fractions down to 10~ 14 are 
displayed for reference purposes. However, for reaction networks coupled to hydro- 
dynamics only mass fractions larger than say ~ 10 2 — 10~ 3 are likely to have sig- 
nificant influence on the hydrodynamics. Thus, the largest discrepancies between 
PE and implicit mass fractions in Fig. [TTT a). and in the other examples that will be 
discussed, imply uncertainties in the total mass being evolved by the network that 
would be irrelevant in a coupled hydrodynamical simulation. 

Timestepping for various integration methods is illustrated in Fig. [TTT b). The 
PE calculation required 3941 total integration steps while the implicit code required 
only 600 steps, but this factor of 6.5 timestepping advantage is offset substantially 
by the expectation that for a 16-isotope network an explicit calculation should be 
about 3 times faster than the implicit calculation for each timestep (Table |3]). Thus 
we conclude that for fully-optimized codes the implicit calculation would be per- 
haps twice as fast for this example. Since the semi-implicit timestepping curve 
was reproduced from another reference the exact number of integration steps is not 
available, but a comparison of the curves in Fig. \TT\b) suggests that an optimized 
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Fig. 11. Calculation for constant T = 5 x 10 9 K (T 9 = 5) and p = 1 x 10 7 gem" 3 alpha 
network (16 isotopes, 48 reactions, and 19 reaction groups). Initial equal mass fractions of 
12 C and 16 0, and rates from REACLIB (191 and Ref. (28j were used, (a) Mass fractions, 
(b) Integration timesteps for implicit code Xnet ll20ll . semi-implicit code YASS |29l (repro- 
duced from from Ref. lfl2V ). asymptotic plus PE (present work), QSS plus PE calculation 
reproduced from Ref. [12], QSS (present work), and asymptotic (present work). 
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Fig. 12. Fraction of reactions that are treated as being in equilibrium as a function of in- 
tegration time for the calculation of Fig. [TT] Also shown is the fraction of isotopes that 
become asymptotic in the PE calculation. 

partial equilibrium code is likely to be at least as fast as the semi- implicit YASS 
code for this example. 

The Asy + PE and implicit timesteps are many orders of magnitude better 
at late times than those of the purely asymptotic (labeled Asy) and purely quasi- 
steady-state (labeled QSS) calculations. The reason is clear from Fig. [121 For times 
later than ~ log? = —6, significant numbers of reaction groups come into equi- 
librium [as determined by the criteria of Eq. (1211)1 . which asymptotic and QSS 
methods alone are not designed to handle. The earlier application by Mott lfT2l of 
a QSS plus PE calculation for this same system is seen to lag behind the current 
implementation of asymptotic plus partial equilibrium in timestepping at late times 
by a factor of 1000 or more. In fact, the timestepping reproduced from Ref. |[T2l for 
QSS + PE is only a little better than that of the pure QSS results from the present 
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Fig. 13. Asymptotic, partial equilibrium, and implicit calculations for an alpha network 
under Type la supernova conditions, (a) Mass fractions, (b) Integration timesteps. (c) The 
hydrodynamic profile, (d) Fraction of isotopes that become asymptotic and fraction of net- 
work reactions that become equilibrated. Solid green indicates the fraction of asymptotic 
isotopes in the PE calculation; dashed green indicates the fraction of asymptotic isotopes 
in the purely asymptotic calculation. 

paper, and clearly is not competitive with that of either the implicit or semi-implicit 
calculations, or the present asymptotic plus partial equilibrium result. Thus we see 
in this example that the huge speed advantage of implicit methods relative to pure 
asymptotic or quasi-steady-state methods in the approach to equilibrium that was 
demonstrated in §18.41 has been essentially erased by a proper treatment of partial 
equilibrium in the explicit integrations. In addition, the present implementation of 
asymptotic plus PE methods is seen to be much faster than previous applications of 
explicit partial equilibrium methods to this network. 



9.2 Example with a Hydrodynamical Profile 



The preceding example employed an alpha network at extreme but constant 
temperature and density. A partial equilibrium calculation using a hydrodynamical 
profile with the dramatic temperature rise characteristic of a burning wave in a Type 
la supernova simulation is illustrated in Fig. [[3] These are challenging conditions 
for the reaction network. During the thermonuclear runaway the temperature in- 
creases by 3.6 billion K in only 2 x 10~ 8 s, corresponding to a rate of temperature 
change 1.7 x 10 17 K/s, and the fastest and slowest rates in the network differ by 
approximately 10 orders of magnitude. Again we see that the partial equilibrium 
timestepping is competitive with that of the implicit code. The partial equilibrium 
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method required 596 total integration timesteps and the implicit code required 553 
steps, so an optimized partial equilibrium calculation would be perhaps three times 
faster than the corresponding implicit calculation because of its speed advantage in 
computing each timestep. At late times the timestepping for both the implicit and 
partial equilibrium calculations is again found to be orders of magnitude larger 
than that from the purely asymptotic calculation. The reason can be seen from 
Fig. fT3T d): at late times the network is very strongly equilibrated and can be in- 
tegrated efficiently with an explicit method only if the equilibrating reactions are 
identified and removed from the numerical integration. 

The synergism of the asymptotic and partial equilibrium approximations is il- 
lustrated by Fig. EE2d), where we see that in the asymptotic plus PE calculation 
isotopes become asymptotic much earlier than in the purely asymptotic calculation. 
Thus the much larger timestep that the PE plus asymptotic calculation takes relative 
to the purely asymptotic calculation in the time interval log? ~ —7 to log? ~ — 3 
is enabled by a combination of reduced stiffness because of the PE approximation 
and reduced stiffness because more isotopes become asymptotic at early times. 
This synergism is common and we find typically that replacing the stiffest parts of 
a network by a complementary set of algebraic constraints exploiting both asymp- 
totic and partial equilibrium conditions is much more effective than either set of 
constraints used alone. 

9.3 Synopsis of Partial Equilibrium Results 

The examples shown above are representative of various problems that have 
been investigated with the new asymptotic plus PE methods. Our general conclu- 
sion is that for alpha networks under the extreme conditions corresponding to a 
Type la supernova explosion, the asymptotic plus partial equilibrium algorithm is 
capable of timestepping that is orders of magnitude better than asymptotic or QSS 
approximations alone when the system approaches equilibrium. These timesteps 
typically lie in the same ballpark as those from current implicit and semi-implicit 
codes, once the partial equilibrium algorithm has removed the fast timescales as- 
sociated with approach to equilibrium. Because the explicit methods can compute 
the timestep more efficiently, we find that in many cases they project to be at least 
as fast as current implicit codes, even for the relatively small networks used as 
examples here. 

10 Improving the Speed of Explicit Codes 

The preceding discussion has assumed that the relative speed of explicit and 
implicit methods can be compared — even if they are presently implemented with 
different levels of optimization — by comparing the ratio of integration timesteps 
required to do a particular problem with each method multiplied by a ratio of times 
to execute a step with each method. The assumption of this comparison is that once 
the matrix overhead associated with implicit methods is factored out, the rest of 
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the work required in a timestep is similar for explicit and implicit methods. This 
is useful for a first estimation, and is sufficient for the primary purpose of this pa- 
per, which is to determine whether explicit methods can take stable timesteps that 
are even in the same ballpark as implicit methods. However, the results presented 
here should not be overinterpreted. Because of different levels of optimization for 
present implicit and explicit timestepping algorithms, and possible differences re- 
sulting from choices of integration tolerance parameters for implicit and explicit 
codes, the present comparisons of integration timesteps are probably uncertain by 
at least an order of magnitude. Thus, the essential message of the present results 
is that, contrary to most previous claims, explicit methods can compete favorably 
with implicit methods for a range of extremely stiff problems. A definitive compar- 
ison of whether implicit or explicit methods are faster for specific problems, and by 
how much, must await more complete development and optimization for the new 
explicit methods. 

As has been noted, the explicit timestepping algorithm employed in this pa- 
per is serviceable but is not very optimized. Thus, it is possible that the number 
of integration steps required by our explicit methods has been overestimated for 
the examples that have been discussed. But in addition to this obvious potential 
for optimization, there are at least two additional reasons why we may not yet be 
realizing the full capability of explicit methods. 

(1) For systems near equilibrium, this paper has emphasized partial equilib- 
rium methods in conjunction with asymptotic approximations. But we have also 
presented evidence that quasi- steady- state (QSS) methods give results comparable 
to asymptotic methods, often with timesteps that can be somewhat larger. Thus, it 
is possible that a QSS plus partial equilibrium approximation could give even bet- 
ter timestepping than the asymptotic plus partial equilibrium examples shown here. 
This possibility has not yet been investigated. 

(2) Potentially of most importance, the observation that algebraically- stabilized 
explicit methods may be viable for large, extremely-stiff networks changes the rules 
of optimization for large reaction networks. For standard implicit methods, the 
most effective optimizations improve the numerical linear algebra, since in large 
networks the bulk of the computing time for implicit methods is spent in matrix 
inversions. But with algebraically- stabilized explicit methods there are no matrix 
inversions and the majority of the time is spent in computing the rates at each step. 
Faster computation of reaction rates has only a small impact on the speed of an 
implicit code for large networks, so there has been little incentive to worry about 
this before. But now we see that our new explicit methods can gain even more in 
speed relative to implicit methods by increasing (through software or hardware) the 
efficiency for computing rates at each timestep. These potential speed gains are sep- 
arate from, and in addition to, those accruing from the absence of matrix inversions 
in explicit methods that have been emphasized in preceding sections. 

Let us cite a simple example illustrating this second point. We have found 
in a single-zone hydrodynamics calculation under Type la supernova conditions 
that using a simple rate-interpolation scheme to avoid recomputing rates for every 
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hydrodynamical timestep permits the speed for explicit asymptotic integration for a 
150-isotope network to be increased by a factor ~ 20 in the region of strong silicon 
burning. 

Also of interest in this context is whether the architecture of modern multicore 
processors with associated GPU accelerators can be exploited to compute rates 
more efficiently for large-scale simulations of reaction networks coupled to fluid 
dynamics. We reiterate that such optimizations could be used for implicit integra- 
tions too, but they will not increase their speeds by as much because rates are a 
smaller part of the computing budget for implicit integration of a large network. 

11 Extension of Partial Equilibrium to Larger Thermonuclear Networks 

Well away from equilibrium, asymptotic and QSS approximations can compete 
with current implicit codes even for extremely stiff systems. There are important 
cases where the system does not equilibrate strongly during the evolution of most 
physical interest (for example, those described in § S8.H - [QT) . Even in problems 
where equilibrium is important overall, there will be many hydrodynamical zones 
and timesteps for which equilibrium plays a small role. In these situations, the 
asymptotic and QSS approximations alone may be adequate to solve large networks 
efficiently. However, to compete systematically with implicit solvers across a broad 
range of problems, asymptotic and QSS approximations must be augmented by PE 
methods to remain viable in those zones where equilibrium is important. 

The examples shown here have demonstrated that asymptotic (and presum- 
ably QSS) plus partial equilibrium methods can solve stiff thermonuclear alpha 
networks with timestepping similar to that for implicit or semi-implicit codes, and 
accuracy more than sufficient to couple such an algorithm to fluid dynamics simu- 
lations. Explicit methods scale linearly and hence more favorable with network size 
than implicit methods, so their relative advantage grows for larger networks. Thus, 
for explicit methods to realize their full potential the previous partial equilibrium 
examples must be extended to include very stiff networks containing hundreds (or 
more) species. That work is ongoing and will be reported in future publications, but 
we note that recent encouraging results in this direction are discussed in Ref. IfTTl . 

12 Summary 

Previous discussions of integrating extremely stiff reaction networks have usu- 
ally concluded that such systems can be integrated only by implicit schemes, be- 
cause explicit schemes are unstable for timesteps large enough to be efficient. Nev- 
ertheless, an alternative to implicit integration for stiff systems is to modify the 
equations to be integrated using approximate algebraic solutions to reduce the stiff- 
ness, and then to integrate the resulting equations numerically by explicit means. 
For example, explicit asymptotic and QSS methods have had some success in in- 
tegrating moderately stiff systems in various chemical kinetics problems Hill 1112.11 . 
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However, these same methods have been found wanting for the extremely stiff 
systems characteristic of astrophysical thermonuclear networks, both because they 
failed to give reliable results, and because the integration timestep was found not 
to be competitive with implicit methods, even if the results had been correct H1I12H . 
We have presented evidence strongly challenging all of this conventional wisdom. 

The key to realizing this new view of explicit integration for stiff systems is 
the understanding that in reaction networks there are (at least) three fundamental 
sources of stiffness: 

(1) Negative populations, which can evolve from an initially small positive popu- 
lation if an explicit timestep is too large. 

(2) Macroscopic equilibration, where the right sides of the differential equations 
dY = F = F + — F in Eq. (Q~|) approach a constant. 

(3) Microscopic equilibration, where the net flux in specific forward-reverse re- 
action pairs (/+ — on the right side of Eq. (OQ) tends to zero. 

These distinctions are crucial because the algebra required to stabilize the explicit 
solution depends on which forms of stiffness are present in a network. The first 
two sources of stiffness are handled well by asymptotic and QSS methods, but 
microscopic equilibration can be handled efficiently using explicit methods only 
by employing partial equilibrium approximations. 

By comparing the timesteps required to integrate problems using explicit and 
implicit methods, and by using an implicit backward-Euler code to estimate the 
relative speed for computing a timestep by explicit and implicit means in networks 
of various sizes, we have shown that algebraically- stabilized explicit integration 
can give correct results and competitive timesteps for even the stiffest of networks. 

(1) For smaller networks without significant microscopic equilibration, we find 
systematic evidence that the explicit methods developed here can give inte- 
gration speeds comparable to that for implicit integration. 

(2) For some cases without much equilibration we find evidence that explicit 
methods might be capable of outperforming implicit ones. We have estimated 
that for equivalently-optimized codes the QSS method used here could be 
faster than the implicit method used here by factors of ~ 10 — 15 for the nova 
simulation of §18.2[ ~ 10 — 15 for the tidal supernova alpha network of § 18.31 
and ~ 5 for the tidal supernova 365-isotope network of ^83] (though we have 
noted that these estimates are uncertain by at least an order of magnitude). 

(3) For those cases where equilibration is significant the asymptotic and QSS 
methods give correct results, but when used alone are much too slow to com- 
pete with implicit methods. However, we presented evidence in §(9] and §QT| 
that if explicit partial equilibrium methods are used to remove the stiffness as- 
sociated with the approach to equilibrium, the explicit methods again exhibit 
speeds that are comparable to or greater than those for implicit methods. 

In addition, we have argued in § fl0l that these estimates may not yet represent the 
best speeds for the new explicit methods. Because 

(1) our timestepping algorithm is not yet optimal, 

(2) we have yet to investigate whether (as might be expected) QSS + partial equi- 
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librium can take larger timesteps than asymptotic + partial equilibrium, and 
(3) an increase (algorithmically or through hardware) in the speed for computing 

rates should preferentially benefit explicit methods relative to implicit ones, 
it is possible that the speed of explicit methods relative to implicit ones might be 
increased substantially relative to the present results. 

Finally, the present results owe a considerable debt to previous work H1I11I12H . 
However, we find that our versions of asymptotic and quasi-steady-state methods 
are much more accurate, and our versions of partial equilibrium methods are much 
faster, than those of previous authors when applied to extremely stiff astrophysical 
thermonuclear networks. The test cases presented here represent a mix of quite 
extreme conditions, with temperature changes as rapid as ~ 10 17 K/s, differences in 
fastest and slowest network timescales as much as 10-20 orders of magnitude, and 
the fraction of equilibrated reactions often changing rapidly between 0% and 100%. 
Hence, we conjecture that these methods may be even more useful for problems 
where the conditions are not quite as extreme as for astrophysical thermonuclear 
environments. 

13 Conclusions 

This paper demonstrates that algebraically- stabilized explicit integration is ca- 
pable of timesteps competitive with those of implicit methods for various extremely- 
stiff reaction networks. Since explicit methods can execute a timestep faster than an 
implicit method in a large network, our results suggest that algebraically- stabilized 
explicit algorithms may be capable of performing as well as, or even substantially 
outperforming, implicit integration in a variety of moderate to extremely stiff ap- 
plications. Because of the linear scaling with reaction network size for explicit 
methods, this fundamentally new view of explicit integration for stiff equations 
is particularly important for applications in fields where more realistic — and there- 
fore larger — reaction networks are required for physical simulations. Arguably, this 
means almost all scientific and technical disciplines, since the sizes of reaction net- 
works being used in simulations to this point have been dictated more often by what 
was feasible than by what was physical. Of particular significance is that these new 
explicit methods might permit coupling of more physically-realistic reaction kinet- 
ics to fluid dynamics simulations in a variety of disciplines. 
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A Appendix: Reaction Group Classification 



Applying the principles discussed §15.4| to the reaction group classes in Table[2]gives 
the following partial equilibrium properties of reaction group classes for astrophys- 
ical thermonuclear networks. 

Reaction Group Class A (a # b) 

Source term: = -k f y a + k r y b Constraints: y a + y b = ci = y a + y^ 
Equation: ^ = by a + c b = -k f c = k T Solution: y a (?) = y° a e ht - § ( 1 - e ht ) 
Equil. solution: y a = — | = j| Equil. timescale: T = j = j- f 

Equil. tests: < £, (i = a,b) Equil. constraint: ^ = j| 
Other variables: y b = c\—y a 

Progress variable: X =y° a -y a y a =y a -X y b =y b ) + X 



Reaction Group Class B(a + bf i c) 

dy a 
dt 



Source term: % = -k f y a y b +k r y c 



Constraints : y b -y a = a =y° h -y° a y b +y c = c 2 = y b +y° 
Equation: ^ = ay a + by a + c a = —kf b = —(c\kf + k b ) c = k r (c 2 — c\] 
Solution: Eq. (fT6l) Equil. solution: Eq. (flTT) Equil. timescale: Eq. (l20l) 
Equil. tests: fc^i < £, (i = a,b, c) Equil. constraint: 2a£ = | 
Other variables: y b = c\ +y a y c = c 2 — y b 

Progress variable: X=y° a -y a y a =y a -X y b =y b -X y c =y° c + X 
Reaction Group Class C (a + b + c d) 

Source term: 4ja = -k f y a y b y c + k r y d Constraints: y a - y b = c\ = y a - y° h 

3 (y a + yb + y c ) +yd = c 3 = 3 (y° +y° b + y° c ) + y° d y a - y c = c 2 =y a ~y Q c 

Equation: -J* = ay 2 a + by a + c 

a= —k f y^ + k i (c\+C2) b = —(k f c\C2 +k T ) c = (c 3 + |ci + ^2)^ 
Solution: Eq. (fT6l) Equil. solution: Eq. (PT7T) Equil. timescale: Eq. (1201) 

Equil. tests: < ej (i = a,b, c,d) Equil. constraint: ™£ = | 
Other variables: y b =y a - C \ y c =y a -c 2 y d = c 3 -y a + \(ci +c 2 ) 
Progress variable: X=y° a -y a y a =y a -X y b =y° b -X y c =y Q c -X 
y d = X+y° d 



c 
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Reaction Group Class D(a + bv^c + d) 

Source term: = -k f y a y h + k r y c y d Constraints: y a - y b = c\ = y a - y Q b 

ya+yc = c2=y° a +y° c y a +yd = c3=y° a +y d 

Equation: -J? = ay 2 a + by a + c a = k r — k f b = —k x {c 2 + C3) + k f c\ c = k r c 2 C3 
Solution: Eq. (fT6l) Equil. solution: Eq. (PT7T) Equil. timescale: Eq. (1201) 
Equil. tests: J^izM < £, (z = a,b, c,d) Equil. constraint: ^ = | 
Other variables: y b = y a - c\ y c = c 2 - y a v rf = c 3 - y a 

Progress variable: X=y a -y a y a = y a ~ & yb=yl~k y c = y° c + A yd = 
y d + X 

Reaction Group Class E(a + bv^c + d + e) 

Source term: ^ = -kfy a y b + k r y c y d y e 

Constraints: y a + ±(y c +y d +y e ) = a = y° a + + y d +y° e ) 

ya-yb = C2=y a -y Q b y c -y d = c^=y° c -y d y c -y e = c 4 =y c -y° e 

Equation: = ay 2 a + by a + c 

a= (3ci -y[j)/c r -/cf b = c 2 k { - (oc/3 + ay4-/3y)/c r c = /c r a/3y 

a = Ci + ^(C3+C4) P = Ci - §C3 + \c$ 7=Ci+^C3-§C4 

Solution: Eq. (fT6l) Equil. solution: Eq. (flTT) Equil. timescale: Eq. (l20l) 
Equil. tests: < £, (z = a,b,c,d,e) Equil. constraint: ^ggj = | 

Other variables: y b = y a - c 2 y c = cx-y a yd=P~y a y e = J-y a 
Progress variable: X=y° a -y a y a =y a -X y b =y° b -X y c =y° c + X 
yd=y° d + k y e = y° + A 

In the equilibrium test condition we have allowed the possibility of a different £, 
for each species i but in practice one would often choose the same small value £ for 
all i. The results presented in this paper have used £, = 0.01 for all species. 
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